library(Seurat)
library(tidyseurat)
library(harmony)
library(SingleR)
library(ggpubr)
library(ggrepel)
library(tidyverse)
source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

sobj <- read_rds('mission/ye_AH-cOME/ah9.rds')

sobj |>
  DotPlot(c('S.Score', 'G2M.Score'))

# annotate with markers in Wang-Ye paper
ye.bcell.type <- list(
  naive.B = c('TCL1A','SELL','FCER2'),
  FCRL4.MBC = c('SOX5','FCRL4','FCRL5','PTPN1','PLAC8','ITGAX'),
  classic.MBC = c('TNFRSF13B','BANK1','CD27'),
  LZ.GC = c('CD83','BCL2A1'),
  DZ.GC = c('CD38','SUGCT','AICDA'),
  cycling = c('HMGB2','MKI67','TUBA1B'),
  plasmablast = c('XBP1','MZB1','JCHAIN'),
  CD21low.MBC = c('COX5A','MYL6','UQCRH','COTL1')
)

sobj |>
  DotPlot(ye.bcell.type, cluster.idents = T, cols = 'RdYlBu',
          group.by = 'seurat_clusters') +
  rotate_x_text()

sobj <-
"
RNA_snn_res.0.8,manual_wang
1,Naive.B
2,DZ.GC
3,cycling
4,cycling
5,FCRL4.MBC
6,plasmablast
7,CD21low
8,LZ.GC
9,cycling
10,DZ.GC
11,classical.MBC
12,plasmablast
" |> read_delim() |>
  mutate(RNA_snn_res.0.8 = as_factor(RNA_snn_res.0.8)) |>
  left_join(sobj, y = _)

sobj |>
  DimPlot(group.by = 'manual_wang', cols = 'Set3', label = T) +
  ggtitle('Adenoid B cells subsets')

sobj |>
  DotPlot(list_c(ye.bcell.type), group.by = 'manual_wang', cluster.idents = T) +
  RotatedAxis()

sobj |>
  mutate(cOME = ifelse(cOME == 1, 'cOME', 'Ctrl')) |>
  write_rds('mission/ye_AH-cOME/ah9.rds')

# annotate with ref data ------
hpca.immu <- read_rds('00_util_scripts/ref/hpca_immune.rds')

sobj %<>% mark_cell_type_singler(hpca.immu, fine_label = T,
                                 new_label = 'hpca.fine')

sobj |> DimPlot(group.by = 'hpca.fine')

sobj <- sobj |> mutate(hpca.fine = case_match(hpca.fine,
                                              'B_cell:Germinal_center' ~ 'GC B',
                                              'B_cell:Memory' ~ 'Memory B',
                                              'B_cell:Naive' ~ 'Naive B',
                                              .default = 'Plasma cell'))

monaco <- celldex::MonacoImmuneData()

sobj %<>% mark_cell_type_singler(monaco, fine_label = T,
                                 new_label = 'mona.fine')

sobj |> DimPlot(group.by = 'mona.fine')

sobj |> FeaturePlot(c('FCRL4','AICDA','CD38'))

latent.come <- sobj |>
  annotate_latents('B.cell', logfc.thres = 1.5)

# compare cell fraction between group -----------
## cOME vs ctrl --------
frac.come <- sobj@meta.data |>
  discov_frac_change(cOME, manual_wang, cOME, Ctrl)

frac.come |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col() +
  scale_fill_manual(values = c('blue','red','grey')) +
  theme_pubr() +
  coord_flip() +
  labs(title = 'B cell subset fraction changes in AH + cOME vs AH patients')

sobj@meta.data |>
  dplyr::count(orig.ident, manual_wang, cOME) |>
  group_by(orig.ident) |>
  reframe(cOME ,manual_wang, fraction = n / sum(n)) |>
  mutate(subtype = manual_wang, .keep = 'unused') |>
  left_join(frac.come[,c(1,5)]) |>
  mutate(subtype = fct_reorder(subtype, log2fc_frac, .desc = T),
         cOME = ifelse(cOME == 'Ctrl', 'AH','AH+cOME') |> fct_relevel('AH')) |>
  ggplot(aes(cOME, fraction, fill = cOME)) +
  stat_summary(geom = 'col') +
  geom_jitter(width = .01, height = 0) +
  facet_wrap(~subtype, scales = 'free_y') +
  theme_pubr() +
  scale_fill_manual(values = c('skyblue','red')) +
  labs(title = 'B cell subset fraction changes by sample',
       fill = 'group', x = 'group')

## RR vs GG+GR --------
frac.rr <- sobj@meta.data |>
  mutate(RRvG = ifelse(`IGHG1-G396R` == 'RR', 'RR', 'other')) |>
  discov_frac_change(RRvG, manual_wang, RR, other)

frac.rr |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col() +
  scale_fill_manual(values = c('blue','red','grey')) +
  theme_pubr() +
  coord_flip() +
  labs(title = 'B cell subset fraction changes in RR vs GG+GR patients')

## IT vs II -----
frac.it <- sobj@meta.data |>
  discov_frac_change(`FCGR2B-I232T`, manual_wang, IT, II)

frac.it |>
  ggplot(aes(subtype, log2fc_frac, fill = type)) +
  geom_col() +
  scale_fill_manual(values = c('blue','red','grey')) +
  theme_pubr() +
  coord_flip() +
  labs(title = 'B cell subset fraction changes in IT vs II patients')

# DEG ---------
## between AH+COME vs AH --------
Idents(sobj) <- 'manual_wang'

Idents(sobj) <- 'hpca.fine'

come.type.list <- unique(Idents(sobj))

come.type.list <- come.type.list |>
  set_names(come.type.list)

deg.subsets <- come.type.list |>
  map(\(x)FindMarkers(sobj, ident.1 = 'cOME', group.by = 'cOME',
                      subset.ident = x) |>
        as_tibble(rownames = 'gene'), .progress = T) |>
  list_rbind(names_to = 'subtype')

sdeg.subsets <- deg.subsets |>
  select(-p_val) |>
  filter(p_val_adj < .05) |>
  write_csv('mission/ye_AH-cOME/deg.subsets.csv')

sdeg.subsets <- sdeg.subsets |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj))

volc.list <- come.type.list |>
  map(\(x)filter(sdeg.subsets, subtype == x, str_starts(gene, 'ENSG|LINC', negate = T)) |>
        plot_bill_volc('AH+cOME') +
        ggtitle(x))
  
volc.list[[8]]

## between G2R --------
### RR vs GG+GR
deg.subsets.g2r <- come.type.list |>
  set_names(come.type.list) |>
  map(\(x)FindMarkers(sobj, ident.1 = 'RR', group.by = 'IGHG1-G396R',
                      subset.ident = x) |>
        as_tibble(rownames = 'gene'), .progress = T) |>
  list_rbind(names_to = 'subtype')

sdeg.subsets.g2r <- deg.subsets.g2r |>
  select(-p_val) |>
  filter(p_val_adj < .05) |>
  write_csv('mission/ye_AH-cOME/deg.subsets.g2r.csv')

sdeg.hpca.rrvgx <- deg.subsets.g2r |>
  select(-p_val) |>
  filter(p_val_adj < .05) |>
  write_csv('mission/ye_AH-cOME/deg.hpca.rrvgx.csv')

### GG vs GR+RR
deg.subsets.gg <- come.type.list |>
  set_names(come.type.list) |>
  map(\(x)FindMarkers(sobj, ident.1 = 'GG', group.by = 'IGHG1-G396R',
                      subset.ident = x) |>
        as_tibble(rownames = 'gene'), .progress = T) |>
  list_rbind(names_to = 'subtype')

sdeg.subsets.gg <- deg.subsets.gg |>
  select(-p_val) |>
  filter(p_val_adj < .05) |>
  write_csv('mission/ye_AH-cOME/deg.subsets.gg.csv')

sdeg.hpca.rxvgg <- deg.subsets.gg |>
  select(-p_val) |>
  filter(p_val_adj < .05) |>
  mutate(avg_log2FC = -avg_log2FC) |>
  write_csv('mission/ye_AH-cOME/deg.hpca.rxvgg.csv')

## RR vs GG
deg.subsets.rrgg <- come.type.list |>
  set_names(come.type.list) |>
  map(\(x)FindMarkers(sobj, ident.1 = 'RR', ident.2 = 'GG', group.by = 'IGHG1-G396R',
                      subset.ident = x) |>
        as_tibble(rownames = 'gene'), .progress = T) |>
  list_rbind(names_to = 'subtype')

sdeg.hpca.rrvgg <- deg.subsets.rrgg |>
  select(-p_val) |>
  filter(p_val_adj < .05) |>
  write_csv('mission/ye_AH-cOME/deg.hpca.rrvgg.csv')

## between I2T --------
deg.subsets.i2t <- come.type.list |>
  set_names(come.type.list) |>
  map(\(x)FindMarkers(sobj, ident.1 = 'IT', group.by = 'FCGR2B-I232T',
                      subset.ident = x) |>
        as_tibble(rownames = 'gene'), .progress = T) |>
  list_rbind(names_to = 'subtype')

sdeg.subsets.i2t <- deg.subsets.i2t |>
  select(-p_val) |>
  filter(p_val_adj < .05) |>
  write_csv('mission/ye_AH-cOME/deg.subsets.i2t.csv')

# try assign IG isotype B cell
ig.isotype <- sdeg.subsets |>
  filter(str_detect(gene, '^IGH[D|M|A|E|G]')) |>
  pull(gene) |> unique()

cb.isotype <- sobj |>
  get_abundance_sc_long(ig.isotype, exclude_zeros = T)

pure.isotype <- cb.isotype |>
  dplyr::count(.cell) |>
  filter(n == 1) |>
  left_join(cb.isotype) |>
  mutate(isotype = .feature, .cell = .cell, .keep = 'none')
  
pure.isotype <- sobj@meta.data |>
  as_tibble(rownames = '.cell') |>
  right_join(pure.isotype)

pure.isotype |>
  mutate(cOME = ifelse(cOME == 'cOME', 'AH+cOME','AH')) |>
  ggplot(aes(cOME, fill = isotype)) +
  geom_bar(position = 'fill') +
  coord_flip() +
  scale_fill_brewer(palette = 'Paired') +
  theme_pubr() +
  labs(title = 'B cell isotype fraction of AH+cOME and AH patients',
       y = 'Fraction',
       x = 'Group')

pure.isotype |>
  ggplot(aes(`IGHG1-G396R`, fill = isotype)) +
  geom_bar(position = 'fill') +
  coord_flip() +
  scale_fill_brewer(palette = 'Paired') +
  theme_pubr() +
  labs(title = 'B cell isotype fraction of AH patients with IGHG1-G396R genotype',
       y = 'Fraction',
       x = 'IGHG1-G396R genotype')

pure.isotype |>
  ggplot(aes(`FCGR2B-I232T`, fill = isotype)) +
  geom_bar(position = 'fill') +
  coord_flip() +
  scale_fill_brewer(palette = 'Paired') +
  theme_pubr() +
  labs(title = 'B cell isotype fraction of AH patients with FCGR2B-I232T genotype',
       y = 'Fraction',
       x = 'FCGR2B-I232T genotype')

# key cytokine violin --------
sobj |> VlnPlot(c('TNF','IL1B','IL4','IL6','CXCL8','IL17A'),
                pt.size = 0)

sobj |>
  filter(str_detect(manual_wang, 'MBC')) |>
  VlnPlot('TNF', group.by = 'IGHG1-G396R') +
  labs(title = 'TNF expression in memory B cells from AH patients',
       fill = 'IGHG1-G396R genotype', x = 'IGHG1-G396R genotype')

sobj |>
  filter(str_detect(manual_wang, 'MBC')) |>
  mutate(cOME = ifelse(cOME == 'cOME','AH+cOME','AH')) |>
  VlnPlot('TNF', group.by = 'cOME') +
  labs(title = 'TNF expression in memory B cells from AH patients',
       x = 'Group')

# allele exlcuded G2R B cell in GR individual -------
bc.g2r.call <- read_csv('mission/ye_AH-cOME/ah-come.g2r.bc.csv')

bc.g2r.call <- bc.g2r.call |>
  mutate(.cell = str_c(pid, '_', .cell))

bc.g2r.call |>
  summarise(n = n(), .by = .cell) |>
  filter(n > 1) |>
  left_join(bc.g2r.call) |>
  slice_max(baseq, by = .cell, with_ties = F) |>
  mutate(g2r.AA = case_match(g2r.base, 'C' ~ 'only G',
                               'T' ~ 'only R', .default = 'Both G&R')) |>
  ggplot(aes(g2r.AA, fill = g2r.AA)) + geom_bar()

al.ex.gr <- bc.g2r.call |>
  filter(str_detect(g2r.base, 'X', negate = T), baseq > 10) |>
  left_join(ah.meta) |>
  mutate(.cell = str_c(id, '_', .cell), g2r.base, mapq, baseq,
         .keep = 'none') |>
  slice_max(baseq, by = .cell, with_ties = F) |>
  left_join(x = sobj, y = _) |>
  filter(!is.na(g2r.base))

al.ex.gr |>
  mutate(g2r.aa = ifelse(g2r.base == 'C', 'G', 'R')) |>
  ggplot(aes(g2r.aa, fill = manual_wang)) +
  geom_bar(position = 'fill') +
  scale_fill_brewer(palette = 'Paired') +
  coord_flip() +
  labs(y = 'Fraction', x = 'AA residue on IGHG1-G396R', fill = 'B cell subtype')

sdeg.al.ex.g2r <- al.ex.gr |>
  FindMarkers(group.by = 'g2r.base', ident.1 = 'T') |>
  as_tibble(rownames = 'gene')

sdeg.al.ex.g2r |>
  filter(p_val_adj < .05)

sdeg.al.ex.g2r |>
  mutate(p_val_adj = p_val) |>
  plot_bill_volc(exp_group = 'R') +
  ylab('-log10(p)')

sobj |>
  bill.violin('TNFRSF11B', `IGHG1-G396R`)

# Fc glyco enzyme ----------
glyco.enz <- read_csv('mission/Ig.Fc.glyco.enyzme.csv')

rrvgg <- read_csv('mission/ye_AH-cOME/deg.subsets.rrgg.csv')

rrvgg |>
  filter(gene %in% glyco.enz$gene) |>
  ggplot(aes(subtype, gene, color = avg_log2FC,
             size = -log10(p_val_adj))) +
  geom_point() +
  theme_pubr(legend = 'right', x.text.angle = 45) +
  scale_color_distiller(palette = 'RdYlBu') +
  labs(title = 'Fc glycolysation enzyme in B cell subsets\nRR vs GG')

rrvgx <- read_csv('mission/ye_AH-cOME/deg.subsets.g2r.csv')

rrvgx |>
  filter(gene %in% glyco.enz$gene) |>
  ggplot(aes(subtype, gene, color = avg_log2FC,
             size = -log10(p_val_adj))) +
  geom_point() +
  theme_pubr(legend = 'right', x.text.angle = 45) +
  scale_color_distiller(palette = 'RdYlBu',limits = c(-1.5,1.5)) +
  labs(title = 'Fc glycolysation enzyme in B cell subsets\nRR vs GG+GR')

ggvrx <- read_csv('mission/ye_AH-cOME/deg.subsets.gg.csv')

ggvrx |>
  filter(gene %in% glyco.enz$gene) |>
  mutate(avg_log2FC = -avg_log2FC) |>
  ggplot(aes(subtype, gene, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  theme_pubr(legend = 'right', x.text.angle = 45) +
  scale_color_distiller(palette = 'RdYlBu') +
  labs(title = 'Fc glycolysation enzyme in B cell subsets\nGR+RR vs GG')

sobj |> write_rds('mission/ye_AH-cOME/ah9.rds')

sobj |>
  DotPlot2d('ST6GAL1', group.x = hpca.fine, group.y = `IGHG1-G396R`) +
  RotatedAxis() +
  labs(title = 'ST6GAL1 expression in human adenoid B cell',
       x = 'B cell subset', y = 'IGHG1-G396R genotype') +
  scale_x_discrete(labels = c('GC B','Memory B','Naive B','Plasma cell'))

## hpca type ------
deg.subsets.rrgg |>
  inner_join(glyco.enz) |>
  ggplot(aes(gene, subtype, color = avg_log2FC,
             size = -log10(p_val_adj))) +
  geom_point() +
  geom_text(aes(label = signif(avg_log2FC, 2)), size = 4,
            position = position_nudge(y = .2), color = 'black') +
  theme_pubr(legend = 'right') +
  facet_wrap(~glycosylation, scales = 'free_x') +
  scale_color_distiller(palette = 'RdYlBu') +
  labs(title = 'IgG Fc glycosylation in adenoid B cell: RR vs GG')

deg.subsets.g2r |>
  inner_join(glyco.enz) |>
  ggplot(aes(gene, subtype, color = avg_log2FC,
             size = -log10(p_val_adj))) +
  geom_point() +
  geom_text(aes(label = signif(avg_log2FC, 2)), size = 4,
            position = position_nudge(y = .2), color = 'black') +
  theme_pubr(legend = 'right') +
  facet_wrap(~glycosylation, scales = 'free_x') +
  scale_color_distiller(palette = 'RdYlBu') +
  labs(title = 'IgG Fc glycosylation in adenoid B cell: RR vs GG+GR')

deg.subsets.gg |>
  mutate(avg_log2FC = -avg_log2FC) |>
  inner_join(glyco.enz) |>
  ggplot(aes(gene, subtype, color = avg_log2FC,
             size = -log10(p_val_adj))) +
  geom_point() +
  geom_text(aes(label = signif(avg_log2FC, 2)), size = 4,
            position = position_nudge(y = .2), color = 'black') +
  theme_pubr(legend = 'right') +
  facet_wrap(~glycosylation, scales = 'free_x') +
  scale_color_distiller(palette = 'RdYlBu') +
  labs(title = 'IgG Fc glycosylation in adenoid B cell: RR+GR vs GG')

sdeg.hpca.rxvgg |>
  inner_join(glyco.enz) |>
  ggplot(aes(gene, subtype, color = avg_log2FC,
             size = -log10(p_val_adj))) +
  geom_point() +
  geom_text(aes(label = signif(avg_log2FC, 2)), size = 4,
            position = position_nudge(y = .2), color = 'black') +
  theme_pubr(legend = 'right') +
  facet_wrap(~glycosylation, scales = 'free_x') +
  scale_color_distiller(palette = 'RdYlBu') +
  labs(title = 'IgG Fc glycosylation in adenoid B cell: RR+GR vs GG')

sobj |> DimPlot(group.by = 'hpca.fine')

## find clusters with high fc RR vs GG --------
Idents(sobj) <- 'seurat_clusters'

leiden.list <- sobj$seurat_clusters |> unique()

leiden.rrvgx <- leiden.list |>
  map(\(x)sobj |> FindMarkers(group.by = 'IGHG1-G396R', ident.1 = 'RR',
                              subset.ident = x) |>
        as_tibble(rownames = 'gene') |>
        mutate(cluster = x), .progress = T) |>
  list_rbind()

leiden.rrvgx |>
  filter(p_val_adj < 0.05) |>
  inner_join(glyco.enz) |>
  ggplot(aes(y = gene, x = cluster, color = avg_log2FC,
             size = -log10(p_val_adj))) +
  geom_point() +
  geom_text(aes(label = signif(avg_log2FC, 2)), size = 4,
            position = position_nudge(y = .2), color = 'black') +
  theme_pubr(legend = 'right') +
  facet_wrap(~glycosylation, scales = 'free_y') +
  scale_color_distiller(palette = 'RdYlBu')

leiden.rrvgg <- leiden.list |>
  map(\(x)sobj |> FindMarkers(group.by = 'IGHG1-G396R', ident.1 = 'RR',
                              subset.ident = x, ident.2 = 'GG') |>
        as_tibble(rownames = 'gene') |>
        mutate(cluster = x), .progress = T) |>
  list_rbind()

leiden.rxvgg <- leiden.list |>
  map(\(x)sobj |> FindMarkers(group.by = 'IGHG1-G396R', ident.1 = 'GG',
                              subset.ident = x) |>
        as_tibble(rownames = 'gene') |>
        mutate(cluster = x, avg_log2FC = -avg_log2FC), .progress = T) |>
  list_rbind()

## among subjects --------
list.enz <- glyco.enz |>
  summarise(data = list(gene), .by = glycosylation) |>
  pull(data, name = glycosylation)
  
sobj |>
  mutate(sample = str_c(`IGHG1-G396R`, '_', orig.ident)) |>
  DotPlot(list.enz, group.by = 'sample', cols = 'RdYlBu') +
  RotatedAxis() +
  labs(title = 'IgG Fc glycosylation in adenoid B cell: by individual',
       x = 'Gene', y = 'Individual')

## among cell types ---------
sobj |>
  DotPlot(list.enz, group.by = 'hpca.fine', cols = 'RdYlBu') +
  RotatedAxis() +
  labs(title = 'IgG Fc glycosylation in adenoid B cell: by B cell type',
       x = 'Gene', y = 'Cell type')

# FAM72A-D ------------
fam72x <- sobj |> rownames() |> str_subset('FAM72') |> str_sort()

sobj |> VlnPlot(fam72x, group.by = 'hpca.fine')

sobj <- sobj |>
  mutate(type.fine = case_when(seurat_clusters %in% c(4,7,10,14) ~ 'DZ GC B',
                               seurat_clusters %in% c(12,17) ~ 'LZ GC B',
                               seurat_clusters == 8 ~ 'Memory B',
                               .default = hpca.fine))

g1 <- sobj |>
  filter(str_detect(type.fine, 'GC')) |>
  DimPlot(group.by = 'type.fine')

g2 <- sobj |>
  filter(str_detect(type.fine, 'GC')) |>
  FeaturePlot(fam72x, cols = c('lightgrey','red'), order = T)

g1 + ggtitle('GC B cells') + g2 + plot_layout(widths = 1:2)

g2

sobj |>
  DotPlot(fam72x, group.by = 'hpca.fine', cols = 'RdYlBu', scale = F) +
  labs(title = 'FAM72A/B/C/D expression in human adenoid B cells',
       x = 'Gene', y = 'Cell type')

sobj |>
  DotPlot(fam72x, group.by = 'hpca.fine', cols = 'RdYlBu', scale = F) +
  labs(title = 'FAM72A/B/C/D expression in human adenoid B cells',
       x = 'Gene', y = 'Cell type')

sobj |> bill.violin(fam72x, group.by = hpca.fine) +
  labs(title = 'FAM72A/B/C/D expression in human adenoid B cells',
       fill = 'Cell type', x = 'Cell type', y = 'Expression level')

sobj |> FindMarkers(features = fam72x,
                    group.by = 'hpca.fine', ident.1 = 'GC B')

sobj |> FeaturePlot(fam72x, cols = c('lightgray', 'red'))

ye.bcell.type <- list(
  naive.B = c('TCL1A','SELL','FCER2'),
  classic.MBC = c('TNFRSF13B','BANK1','CD27'),
  LZ.GC = c('CD83','BCL2A1'),
  DZ.GC = c('CD38','SUGCT','AICDA'),
  plasmablast = c('XBP1','MZB1','JCHAIN'),
  S.phase = c('TYMS','PCNA','FEN1'),
  G2.M.phase = c('CDC20','CCNB1','PLK1'))

gc.bcell.type <- list(
  LZ.GC = c('CD83','BCL2A1'),
  DZ.GC = c('CD38','SUGCT','AICDA'),
  S.phase = c('TYMS','PCNA','FEN1'),
  G2.M.phase = c('CDC20','CCNB1','PLK1'))

sobj |>
  filter(str_detect(type.fine, 'GC')) |>
  DotPlot(c(gc.bcell.type, fam72x), cols = 'RdYlBu',
          group.by = 'seurat_clusters') +
  RotatedAxis() +
  scale_y_discrete(labels = 1:6) +
  labs(title = 'GC marker gene & FAM72s expression in adenoid GC B cells',
       x = 'Gene', y = 'Cluster')

## define cell type with cc ----------
sobj <- sobj |>
  mutate(cc.type = case_when(seurat_clusters %in% c(13, 8) ~ 'S phase Memory B',
                             seurat_clusters %in% c(14) ~ 'G2/M phase DZ GC B',
                             type.fine == 'DZ GC B' ~ 'other DZ GC B',
                             .default = type.fine))

sobj |> DotPlot(fam72x, group.by = 'seurat_clusters')

fam72hi.deg <- sobj |>
  FindMarkers(group.by = 'seurat_clusters', ident.1 = 14) |>
  as_tibble(rownames = 'gene')

library(clusterProfiler)

fam72hi.gogse <- fam72hi.deg |>
  filter(p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL')

fam72hi.gogse <- fam72hi.gogse |>
  simplify()

fam72hi.gogse@result |>
  as_tibble() |>
  filter(NES > 0, p.adjust < .05) |>
  slice_sample(n = 10) |>
  mutate(Description = str_wrap(Description, 40) |>
           fct_reorder(NES)) |>
  ggplot(aes(NES, Description, fill = -log10(p.adjust))) +
  geom_col() +
  theme_bw() +
  labs_pubr() +
  scale_fill_distiller(palette = 'Reds', direction = 1) +
  labs(title = 'GO GSEA upregulated pathway in FAM72-high GC B cell') +
  theme(plot.title.position = 'plot')

theme_jpub

## by individual ----------
sobj |>
  DotPlot2d(fam72x, orig.ident, hpca.fine) |>
  pluck('data') |>
  ggplot(aes(group.y, avg.exp, color = group.y)) +
  geom_boxplot() +
  geom_jitter(width = .1, height = 0) +
  facet_wrap(vars(features.plot)) +
  theme_pubr() +
  labs(title = 'FAM72A/B/C/D expression in human adenoid B cells (n=9)',
       x = 'Cell type', y = 'Average expression')

sobj |>
  DotPlot2d(fam72x, orig.ident, cc.type) |>
  pluck('data') |>
  mutate(group.y = fct_reorder(group.y, avg.exp, .desc = F)) |>
  ggplot(aes(group.y, avg.exp, color = group.y)) +
  geom_boxplot() +
  geom_jitter(width = .1, height = 0) +
  facet_wrap(vars(features.plot)) +
  theme_pubr(x.text.angle = 45, legend = 'none') +
  labs(title = 'FAM72A/B/C/D expression in human adenoid B cells (n=9)',
       x = 'Cell type', y = 'Average expression')

# TRPM7 expression with Ab response ------
igh.gene <- rownames(sobj) |> str_subset('^IGH.{1,2}$') |>
  str_subset('P$', negate = T) |>
  str_sort()

sobj |> DotPlot(c('TRPM7', igh.gene), group.by = 'orig.ident',
                cols = 'RdYlBu', cluster.idents = T)

m7.igh.expr <- last_plot() |>
  pluck('data')

m7.igh.expr |>
  filter(features.plot == 'TRPM7') |>
  mutate(id = fct_reorder(id, avg.exp)) |>
  ggplot(aes(id, avg.exp.scaled)) +
  geom_col()

sobj <- sobj |>
  mutate(m7.group = ifelse(orig.ident %in% c('M1025','L0609','C1018','W1101'),
                           'm7-low', 'm7-high'))

m7.hvl <- sobj |>
  FindMarkers(group.by = 'm7.group', ident.1 = 'm7-high') |>
  as_tibble(rownames = 'gene')

m7.hvl |>
  filter(p_val_adj < .05)

m7.hvl |>
  filter(gene %in% c('TRPM7', igh.gene)) |>
  ggplot(aes(gene, avg_log2FC, fill = p_val_adj)) +
  geom_col() +
  theme_bw() +
  labs_pubr() +
  labs(title = 'Expression of IGH genes in adenoid B cells\nTRPM7-high vs TRPM-low')

m7.hvl |>
  filter(gene %in% c('TRPM7', igh.gene)) |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj)) |>
  ggplot(aes(gene, y = 'Average log2FC',
             fill = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point(shape = 21) +
  theme_bw() +
  labs_pubr() +
  scale_size() +
  scale_fill_gradient2(low = 'blue', high = 'red') +
  labs(title = 'Expression of IGH genes in adenoid B cells\nTRPM7-high vs TRPM-low',
       y = NULL)

sobj$mona.fine |> unique()

sobj |>
  DotPlot2d('TRPM7', orig.ident, hpca.fine) +
  labs(x = 'Sample', y = 'B cell subtype', title = 'TRPM7 in adenoid B cells')

sobj@meta.data |>
  distinct(orig.ident, sex, m7.group, `IGHG1-G396R`, `FCGR2B-I232T`, cOME) |>
  arrange(m7.group)
